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,  Introduction  and  Problem  Statement 

Hie  basic  problem  addressed  In  this  paper  is  that  of  providing  an 
accurate  real  time  estimate  of  a  heading  direction  on  a  reference  test 
pad,  Ihis  situation  is  synbolically  depicted  in  Figure  1.  In  this  figure 


Note:  The  main  concepts  in  this  report  were  presented  at  the  1973  AIAA 
Guidance  and  Control  Conference,  Key  Biscayne,  Florida,  20-22  August. 
This  report  is  a  revised  and  corrected  form  of  AIAA  Paper  No.  73-841. 


the  three  dimensional  slab  represents  the  test  pad  to  which  is  fixed  an 
azimuth  measuring  device.  The  function  of  this  device  is  to  measure  the 
angle  between  a  fixed  inertial  direction  and  some  arbitrary  pad  reference. 
Let  us  assume  that  the  quantity  a  is  the  desired  fixed  heading  vector 
and  that  it  is  located  at  an  angle  A  from  the  pad  reference  heading 
vector  o.  Because  the  instrument  is  sensitive  to  other  variables,  we 
cannot  measure  A  directly  but  must  obtain  an  estimate  of  the  a  direc¬ 
tion  by  measuring  the  angle  m  formed  between  the  reference  o  and  the 
indicated  heading  vector  b . 

Much  work  has  been  done  in  modeling  the  error  sources  involved  in 
this  process  [1,2].  It  is  not  the  aim  of  this  paper  to  derive  a  new 
error  model  but  rather  to  demonstrate  a  method  by  which  a  reference 
heading  can  be  obtained  from  raw  gyrocompass  data.  Only  a  sinplified 
error  model  will  be  presented  but  the  proposed  technique  can  be  used 
equally  well  with  any  error  model. 

In  order  to  demonstrate  the  proposed  method  let  us  assume  that  only 
the  platform  North-South  tilt,  T^,  and  East -West  tilt  rate,  T^,  are 
modelled  as  error  sources .  Further  it  is  assumed  that  the  difference 
between  the  angle  m  and  the  angle  A  is  a  linear  function  of  the 
modelled  error  sources.  Thus  we  can  approximate  the  mathematical 
relationship  between  A  and  m  by 

m  -  A  ♦  C;  ♦  C2  .  CD 

Note  that  in  Equation  (1)  the  coefficients  Cj  and  C2  must  be 
determined  along  with  the  heading  A. 
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A  problem  arises  in  that  even  in  the  best  of  error  models  other 
variables  not  modelled  may  affect  the  difference  between  the  measurement 
m  and  the  heading  angle  A.  That  is  to  say,  an  equation  in  the  form  of 
(1)  is  most  likely  not  an  accurate  model  of  the  measurement  process.  In 
order  to  account  for  these  model  inaccuracies  we  can  as sune  that  the 
quantity  A  is  variable  with  time.  If  A  were  assumed  to  be  a  constant 
quantity  along  with  the  coefficients  C j  and  C2  then  there  would  be 
little  hope  of  obtaining  long  term,  real  time  estimates  of  the  heading 
direction.  At  best  all  that  could  be  obtained  would  be  average  estimates 
of  A  for  discrete  time  intervals.  Further,  if  a  standard  Kalman  filter 
were  used  to  try  to  track  these  parameters  the  filter  output  would 
actually  diverge  [3].  Thus  the  problem  is  new  to  estimate  the  constant 

,  where  new 

represents  the  value  of  A  at  the  kth  time  instant. 

A  further  problem  occurs  in  that  the  measurement  of  m  will  also 
br  corrupted  with  measurement  noise.  The  processing  technique  developed 
must  also  account  for  this  fact.  Thus  Equation  (1)  can  be  more 
accurately  represented  by 


error  model  coefficients  and  to  track  the  variable 


+  C,  T, 


NS, 


C2  f, 


EW, 


+  v. 


(2) 


In  the  above  equation  represents  the  noise  and  the  subscript  k 
indicates  that  the  quantities  so  subscripted  are  time  varying  and  that 
their  value  at  the  kth  measurement  is  (*)j.* 
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Our  task  is  thus  to  develop  a  method  by  which  the  measurements 


can  be  processed  so  as  to  yield  accurate  estimates 


of  the  coefficients  Ci  and  C2  while  also  tracking  a  time  variable 
quantity  A^.  We  will  now  proceed  to  derive  a  method  which  can  perform 
the  above  task. 


Method  of  Approach 

Before  proceeding  to  develop  the  method  used  to  solve  the  proposed 
problem,  Equation  (2)  will  be  rewritten  here  in  a  more  general  form. 

To  do  this  the  parameters  to  be  estimated  will  be  written  as  the  column 
vector  x^,  namely 


[\»  ,  C2] 


C3) 


Also  the  measurement  matrix  will  be  defined  as  a  row  vector  M(k) , 
that  is, 

M(k)  -  £l. 

Now  using  the  symbol  y^  to  represent  the  measurement 
Equation  (2)  can  be  written  as 


NY 


EW, 


(4) 


•k  '  M<«  Xk  +  Vk  ‘ 


(5) 
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Equation  (5)  will  be  referred  to  as  the  measurement  equation.  In  the 
above  equation,  vk  is  the  random  measurement  noise  with  zero  mean  and 
covariance  bJv*  vjj  -  R.  In  the  above  It  is  assumed  that  the  neasure- 
mant  equation  is  scalar.  If  more  than  one  heading  measurement  is  to  be 
processed  at  one  time  this  approach  is  still  applicable  though  appro¬ 
priate  changes  must  be  made  to  the  equations  throughout  this  paper. 

It  should  be  noted  here  that  any  error  model  equation  can  be  used 
as  long  as  the  resultant  measurement  equation  is  of  the  same  form  as 
Equation  (5). 

In  general  the  value  of  the  state  vector  at  time  tk+1  can  be 
related  to  its  value  at  time  tk  by  means  of  the  state  transition 
matrix  4>(k+l/k) .  This  allows  us  to  write 

xk+l  "  Hk+l/k)  xk  +  r  Ck)  wk  .  (6) 

Equation  (6)  will  be  referred  to  as  the  state  equation.  In  this  equation 

wk  is  the  system  noise  and  r(k)  is  the  system  noise  coefficient  matrix 

at  time  tk.  The  system  noise  coefficient  matrix  relates  the  effect  that 

the  system  noise  has  upon  the  states.  The  quantity  wk  is  a  random 

variable  with  mean  zero  and  covariance  given  by  Ejwk  wk|  -  q. 

For  the  present  problem  it  is  assumed  that  the  states  are  constant 

between  measurements,  thus  4>(k+l/k)  ■  I.  Because  it  was  assumed  that 

the  system  noise  is  only  in  the  first  state,  r(k)  can  be  given  as 
T 

r  (k)  -  [1,  0,  0].  Thus  all  quantities  in  the  state  and  measurement 
equations  are  specified. 
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The  object  of  the  approach  presented  below  is  to  allow  for  varia¬ 
tions  in  the  parameter  A  by  estimating  the  system  noise  covariance  q. 
This  will  be  done  by  using  the  residual  between  the  actual  measurement 
and  a  predicted  measurement  that  is  based  upon  past  values  of  the 
states  [4,5], 

Assume  that  we  have  processed  the  measurement  yk  and  now  wish  to 
bring  the  filter  and  states  forward  to  time  t^+1>  In  the  processes  of 
going  forward  in  time  we  wish  to  use  that  q  value  which  yields  the 
predicted  residual,  r(k+l/k) ,  that  is  the  most  probable.  In  other 
words,  find  q  according  to 

max  p[r(k+l/k)]  (7) 

q>0 

where  p[  ]  is  the  probability  density  function.  The  q  that 
maximizes  (7)  will  be  denoted  by  q^. 

Let  us  now  define  the  predicted  residual  by 

r(k+l/k)  -  yk+1  -  M(k+1)  x(k+l/k)  .  (8) 

In  this  equation  x(k+l/k)  is  the  estimate  of  the  state  at  time  tk+1 
given  the  measurements  up  to  yk.  The  mean  of  the  predicted  residual 
can  now  be  given  by 
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E[r(k+l/k) ]  -  E[yk+1  -  M(k+1)  x(k*l/k)]  , 

-  M(k+1)  E[xk+1  -  x(k+l/k)]  +  E[vk+1]  , 

-  0  .  (9) 

This  is  by  virtue  of  the  fact  that  the  expected  value  of  the  measurement 
noise  is  assured  to  be  zero  as  is  the  expected  value  of  [xk+1  -  x(k+l/k)], 
the  error  between  the  predicted  state  and  the  true  state. 

Further  the  variance  of  the  residual  can  be  given  by 

E[r(k+l/k)rT(k+l/k)j  - 

Ej[yk+i  “  M(k+l)x(k+l/k)j  [yk+1  -  M(k4l)x(k+1)]TJ  , 
or 

E  |r(k+l/k)rT(k+l/k)  j  - 

E J jMCk+X)  [xk+1  -  x(k+l/k)j  +  vk+1j  [M(k+i)[xk+1  -  x(k  l/k)J  +  vk+1JT|  . 

(10) 

Letting  ek+^  be  the  error  in  the  estimate  at  time  tk+p  that  is, 

Gk+1  *  xk+l  "  *(k+l/k),  equation  (1)  can  be  simplified  to  yield 

E  £r(k+l/k)rT(k+l/k)j  - 
Bj[Mtk*l)cktl  *  VJ  [M(M)ew  *  Vk+1]TJ 
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or 

E[r(k*l/k)rT(k»l/k)l  - 

M(lc»l)  E|cktl  MT(k«)  ♦  E [vk+1  vJjJ 

*  MCk*l)  e|cw  vJjI  *  E|vkn  cjjj  MT(k+l)  . 

The  residual’s  variance  can  now  be  written  as 

E [r (k+l/k) rT(k+l/k) j  - 

M(k+1)  P(k+l/k)  MT(k+l)  +  R(k+1)  ,  (11) 

where  the  state  error  covariance  matrix  is  denoted  is 

P(k+l/k)  -  Ejek+1  ejj)  .  (12) 

Recall  from  the  equations  for  tlie  basic  Kalman  filter  that  the 
system  noise  covariance  enters  into  the  calculation  of  P(k+l/k)  [6], 
This  relationship  can  be  given  by 

P(k+l/k)  -  ♦(k+l/k)  P(k/k)  *T(k+l/k)  +  q  r(k)  rT(k)  ,  (13) 

or  for  this  problem 

/  q  0  o  \ 

P(k+l/k)  -  P(k/k)  +  (  0  0  0  j  .  (14) 

\  0  0  0  / 
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Now  if  we  assume  a  gaussian  distribution  of  the  residuals  for  the 


scalar  measurement  case 


p|r(k+l/k} J 


1  r_ 

'2  02 


(15) 


whei'e 


=  E[r(k+l/k)  rT(k+l/k)J 


(16) 


and 


r2  -  r1 (k+l/k)  r(k+l/k)  .  (17) 

Thus  the  variance  a2  is  a  function  of  q_  and  hence  (15)  is  also 
dependent  on  the  quantity  q. 

This  fact  allows  us  to  minimize  Fquation  (15)  with  respect  to  q 

2  9  fey 2 1 

by  differentiating  it  with  respect  to  o  because  -- — f  is  constant, 

3q 

Proceeding  with  the  differentiation  and  equating  the  result  to  zero 
we  get 


9p(r)  = 


3  o' 


,.2 


1  r2 

7  ^ 


2°  /  2a 2  /2it~c7 


-  e 


1  rz 

7  o7 
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This  implies  that 


I 


which  says  that  the  maximum  probability  occurs  when 


or  equally  when 


r2(k+l/k)  =  E[r?(k+l/k)l  .  (18) 

Now  using  Equation  (13)  in  (11)  we  can  get 

c["r2(k+l/k)]  -  j 

M(k+1)  Hk+1/k)  P  (k/k)  ;T(k+l/k)  MT(k+l) 

+  q  M(k-f-l)  r (k)  rT(k)  MT(k+l)  +  R(k*l)  .  (19) 

Letting 

E[r2(k+l/k)q  «  o]  = 

M(k+1)  <j>(k+l/k)  P(k/k)  ^T(k+l/k)  MT(k+l)  +  R(k+1)  , 


we  can  get  from  (18)  and  (19) 
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if  >  0 


f r2(k+l/k)  -  E[r2(k+l/k)  q  -  0] 
M(k+l)  r(k)  rT(k)  MT(k+l) 


<*k 


(20) 


otherwise . 


where  must  be  >  0  from  the  physical  meaning  of  the  variance. 

Equation  (20)  can  now  be  simplified  for  the  specific  pro-Mem 
described.  This  simplification  allows  us  to  write 


M(k+1)  r(k)  rT(k)  MT(k+l)  - 


[1,tns^ew]i 


[10  0]  |1 


*NS 

L.Vjk 


and 


So 


E[r(k+l/k)q  ■  o]  -  M(k+1)  P(k/k)  MT(k+l)  +  R(k+1) . 


\ 


r2 (k+l/k) -M(k+l)P(k/k)MT(k+l) -R(k+1)  if  >  0 


otherwise. 


(21) 


If  preferred  Equation  (22)  can  also  be  expressed  in  terns  of  the 
measurement  as 
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[Vkf i-IIO^UxOBf l/w]  -R(k+1)  if  >  0 

%  - 

0  otherwise. 

(22) 

Referring  to  Equation  (20)  we  note  that  qk  is  directly  proportional 
to  the  excess  of  the  residual  squared  over  the  predicted  value  of  the 
residual  squared  where  the  predicted  value  is  based  upon  the  assumption 
that  the  system  noise  is  zero.  If  the  residual  squared  is  greater  than 
the  value  predicted  under  the  assumption  of  no  system  noise,  then  a 
positive  value  of  is  generated  which  indicates  that  the  no  noise 
assumption  was  most  likely  false.  This  value  of  q^  is  then  used  in 
later  processings  to  produce  more  consistency  between  the  residuals  and 
their  predicted  statistics. 

The  filter  adapts  to  the  system  noise  level  as  follows.  If 

A 

Equation  (20)  yields  a  non-positive  value  for  qk  then  the  residuals 
are  within  their  la  values  and  the  assumption  of  no  input  noise  was 
probably  true.  The  residuals  are  thus  behaving  according  to  their 
statistics,  are  relatively  small  and  the  filter  is  operating  satis¬ 
factorily.  On  the  other  hand  if  the  residuals  are  larger  than  the  pre¬ 
dicted  la  values ,  the  filter  is  actually  diverging.  This  then 
generates  a  positive  q^  value  which  then  causes  P(k+l/k)  to  increase 
as  seen  from  Equation  (14). 

The  Kalman  gain  can  be  given  by  [7] 

K(k+1)  -  P(k+l/k)  MT(k-0.)  [M(k+1)  P(k+l/k)  MT(k+l)  +  RCk+l)]'1  . 

(23) 
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We  note  that  increasing  P(k+l/k)  causes  the  gain  to  increase.  The 
increased  gain  then  causes  the  filter  to  become  more  sensitive  to  the 
latest  data.  This  allows  the  filter  to  follow  parameter  changes  that 
become  evident  as  later  data  is  processed.  This  is  in  contrast  to  the 
standard  filter  that  tries  to  fit  constant  parameter  values  to  all  of 
the  data  and  which  in  truth  is  biased  towards  earlier  data.  This  u. as 
comes  about  from  the  fact  that  filter  gain  decreases  as  more  and  more 
data  are  prooessed  and  P(k+l/k)  decreases  [5]. 

Die  approach  presented  can  be  used  processing  one  residual  at  a 
time  or  processing  many  residuals  at  once.  This  later  approach  is  more 
statistically  significant  although  does  complicate  the  filter  computa¬ 
tions.  Because  the  value  of  q^  responds  to  large  measurement  noises 
as  well  as  to  large  system  noises ,  a  statistical  approach  to  calculate 
qk  may  be  preferred.  To  process  many  residuals,  we  use  the  sanple 
mean  given  by 


Diis  value  is  then  used  in  the  equations  for  qk  with  appropriate 
modifications. 

Die  main  advantage  of  using  an  adaptive  filter  is  that  it  allows 
the  estimates  of  the  states  to  follow  the  variations  in  the  true  states. 


Thus  the  filter  should  yield  varying  parameters  that  track  the  true 
parameters  accurately,  rapidly  and  with  minimal  computational  burden. 

Figure  2  shows  the  basic  flew  of  information  within  the  filter. 
It  should  be  noted  that  this  differs  from  the  standard  Kalman  filter 


SET 
k  -  1 


INITIAL 

x(k/k) 

P(kA) 


COMPUTE 

x (k+l/k) 
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only  in  the  blocks  that  calculate  and  that  update  P(k+l/k)  if  a 
positive  qN  is  found. 

To  test  the  ability  of  the  adaptive  filter  to  follow  the  variable 
heading  a  sanple  program  was  simulated  using  the  Fortran  IV 
language  on  a  Burroughs  B6 700  computer.  For  this  simulation  it  was 

,i  t  '  ■  ■  . 

assumed  that  the  variation  in  A^  was  sinusoidal  with  a  mean  of  1S.0 
and  a  peak  to  peak  variation  of  6.0.  The  period  of  this  signal  was 
taken  to  be  24  hours .  The  tilt  and  tilt  rate  used  had  a  period  of  1 
hour  and  varied  as  the  sine  and  cosine  respectively.  Both  were  assumed 
to  have  a  zero  mean  with  peak  to  peak  variations  in  the  tilt  of  2.6  and 
in  the  tilt  rate  of  1.6.  The  true  values  of  C!  and  C2  were  chosen 
to  be  0.4  and  0.2. 

For  this  simulation  three  residuals  were  averaged  and  three  measure 
ments  were  processed  simultaneously.  Each  time  a  new  measurement  was 
available  the  oldest  measurement  was  discarded  and  the  three  newest 
measurements  were  processed.  In  this  manner  a  new  estimate  of  A^ 
could  be  obtained  each  time  a  data  point  was  available.  It  was  also 
assumed  that  the  measured  heading,  tilt  and  tilt  rate  were  sampled  every 
10.5  minutes. 

The  results  obtained  when  the  above  case  was  simulated  with  measure 
ment  noise  free  data  being  processed  through  the  filter  are  shown  in 
Figure  3.  For  the  simulation  shown  the  initial  values  of  A^,  Cj  and 
C2  were  chosen  as  12.0,  0.1,  and  0.0.  In  this  figure  we  see  plotted 
the  noise  free  measurement  m^,  the  true  heading  Append  the  heading 
estimate  obtained  from  the  adaptive  filter,  Ap.  Figure  4  shows  the 
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Figure  3.  Results  of  Processing  True  Data  through 
the  Adaptive  Filter 
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Figure  4.  Error  in  Parameter  Estimates  using  True  Data 


error  in  estimating  A^,  C*  and  C2  as  a  function  of  the  nimfcer  of 
points  processed.  As  seen  from  these  two  figures  the  adaptive  filter 
is  able  to  give  quite  accurate  results. 

Ch  the  other  hand  when  the  same  situation  is  simulated  using 
measurements  that  are  corrupted  with  random  noise  the  filter  does  not 
perform  as  well.  In  Figure  5  we  see  the  results  when  the  measurement 
noise  has  a  standard  deviation  of  1.0.  Here  the  noisy  measurement 
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Figure  5.  Results  of  Process ing  Noisy  Data  through 
the  Adaptive  Filter 


is  shown  along  with  the  true  and  estimated  headings  A^  and  Ap.  As 
seen  in  this  figure  the  adaptive  filter  does  not  perform  satisfactorily 


for  this  large  •  measurement  noise.  Even  using  the  averaged  residual 
the  filter  still  reacts  to  the  measurement  noise  as  well  as  to  the 
actual  variations  in  the  heading. 

In  order  for  this  filter  to  behave  properly,  the  noise  level  of  the 
inpdt  signal  must  be  relatively  small.  The  next  portion  of  this  paper 
shows  how  this  secondary  task  can  be  acconplished. 

Noting  that  in  general  the  measurement  signal  will  be  of  a  lower 
frequency  than  the  measurement  noise,  it  appears  that  a  low-pass  filter 
is  needed.  This  type  of  filter  will  allow  low  frequency  signals  to 
pass  almost  wattenuated  while  attenuating  the  higher  frequency  compo¬ 
nents.  Msny  different  types  of  digital  filters  have  been  derived  that 
are  able  to  perfoim  this  task.  A  large  nunber  of  these  are  the  exten¬ 
sion  of  analog  filters  into  the  digital  domain  [8,9]. 

Bach  of  the  available  filter  types  has  advantages  and  disadvantages. 
For  this  application  a  Butterworth  filter  will  be  chosen  because  it  is 
monotonic  in  both  passband  and  in  stopband.  This  type  of  filter  is  able 
to  be  represented  in  digital  form  by  the  following  squared  magnitude 
function 


.  - -a -  ,  (25) 

1  +  tanzn  (uiT/2) 

tan2n  (wcT/2) 

In  this  equation  n  is  the  nuifcer  of  poles  of  the  filter  and  T  is  the 
time  between  sanples.  The  cutoff  frequency,  u  ,  is  the  frequency  at 
which  the  filter  gain  falls  off  to  3  db. 
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Figure  6  shows  the  gain  of  the  digital  Butterworth  filter  as  a 
function  of  frequency  for  various  values  of  n.  While  at  first  it  may 
appear  that  a  large  value  of  n  is  desirable  it  should  be  pointed  out 
that  the  phase  difference  between  the  signal  and  the  filter  output 
increases  with  increasing  n.  Thus  the  value  of  n  should  not  be  any 
larger  than  necessary  if  real  time  data  processing  is  desired. 


Figure  6.  Digital  Butterworth  Filter  Gain  versus  Frequency 
for  a  Variety  of  Poles 

In  order  to  determine  the  minimum  number  of  poles  to  be  used  in  the 
filter  the  desired  operating  criteria  have  to  be  specified.  We  will 
assume  that  the  highest  signal  we  want  to  pass  through  the  filter  is 
1  cycle/hr.,  corresponding  to  the  tilt  and  tilt  rate  effects.  Thus  the 
cutoff  frequency  will  be  1  cycle/hr.  or  27.78  x  10' 5  Hi.  Assuming  that 
that  we  have  a  sanple  every  10.5  minutes  gives  a  sampling  rate  of 
158.73  x  10” 5  Hz .  Further  it  is  specified  that  at  least  a  40  db  drop 


is  required  at  2  cycles /hr.  This  implies  that  the  squared  magnitude 
function  equal  10'4  at  55. 56  x  10’ 5  Hz.  Using  these  numbers  In  Equation 
(25),  n  is  readily  calculated  to  be  3.95 ,  which  implies  that  a  4  pole 
Butterworth  digital  filter  can  perform  the  required  task. 

Proceeding  with  the  design  of  the  filter  we  first  express  the 
general  fourth  order  Eutterworth  filter  as  [10] 

H(s)  - - SSSSlSSS -  .  (26) 

s4  ♦  2.6  s3  ♦  3.4  s2  +  2.6  s  ♦  1 

Replacing  s  by  s/a,  where  a  will  equal  tan  (wcT/2),  Equation  (26) 
can  be  rewritten 


s4  +  2.6  as3  +  3.4  a2s2  ♦  2.6  a3s  +  a4 

In  this  form  the  constant  has  been  adjusted  so  that  the  gain  at  s  ■  0 
is  equal  to  unity. 

Equations  (26)  and  (27)  arc  still  in  analog  form  and  must  now  be 
transformed  into  the  digital  domain.  This  is  accomplished  by  letting 
s  ■  (z-l)/(z+l).  This  transformation  allows  the  filter  transfer 
function  to  be  written  as 

_  a4(z4  +  4z3  +  6z2  ♦  4z  ♦  1) 


az4  +  8z3  +  yZ2  +  6z  ♦  c 


V 


4  ¥ 


The  Greek  symbols  used  in  Equation  (28)  are  defined  by 

a  -  1  +  2.6  a  +  3.4  a2  ♦  2.6  a3  +  a4 
$  ■  -4-5.2a+5.2a3  +  4a4 

Y  -  6  -  6.8  a2  ♦  6  a1*  (29) 

6  -  -  4  +  5.2  a  -  5.2  a3  +  4  a4 
e  •  1  «  2.6  a  +  3.4  a2  -  2.6  a3  ♦  a4 


and  were  obtained  after  algebraic  sinplification  of  the  resulting 
equation  when  the  s  to  z  transform  was  used  in  Equation  (27). 

In  order  to  express  Equation  (28)  in  terns  of  measurement  sample 
times  a  few  more  substitutions  are  required.  First  Equation  (28)  is 
multiplied  by  z‘4/z“4  to  obtain 


a4(l  +  4z-1  +  6z'2  ♦  4z“3  +  z~4' 
a  +  0Z" 1  +  yZ”z  +  6z” 3  +  ez'4 


Recalling  that  the  transfer  function  relates  the  output  of  a  system  to 
its  input,  Equation  (30)  can  be  written  as 

Yf(z) 

Ym« 

where  Y^(z)  is  the  filter  output  and  YJ(l(z)  is  the  filter  input. 


a4(l  +  4z-1  *  6z'z  +  4z'3  +  z~4) 
a  +  ez*1  +  YZ-Z  +  <s  Z" 3  +  eZ'4 
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Equation  (31)  now  relates  the  output  of  the  filter  to  the  input  in  the 
z  plane.  Hie  transformation  to  the  time  stapled  domain  is  now  rather 
straightforward. 

This  process  is  begun  by  expressing  Equation  (31)  in  block  diagram 
form.  Figure  7  is  the  result  of  performing  a  direct  transform  from  the 


I - - — * - Yf  (nT) 

Figure  7.  Block  Diagram  of  Digital  Filter 
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equation  of  the  transfer  function  to  its  block  diagram  form.  Also  in 
this  figure  it  is  shewn  hew  the  z_l  block  relates  to  the  sanples , 
that  is,  a  block  whose  transfer  function  is  z"1  represents  a  delay  of 
one  time  sample.  This  block  diagram  form  now  allows  the  equation  that 
relates  the  filtered  output  to  the  input  to  be  written  as 

«Yf(nT)  +  @Yf(nT-T)  +  yYf(nT-2T)  +  6Yf(nT-3T)  +  eYf(nT-4T)  - 
AYm(nT)  *  BYm(nT-T)  +  CYm(nT-2T)  +  DYm(nT-3T)  +  EYJriIMT)  , 

(32) 

where  the  symbols  used  are  defined  below. 

A  ■  a4 
B  -  4  a4 
C  -  6  a4 
D  ■  4  a4 
E  •  a4 

a  *  1  +  2.6  a  +  3.4  a2  +  2,6  a3  +  a4 

8  -  -  4  -  5.2  a  +  5.2  a3  +  4  a4 

Y  ■  6-6. 8  a2  +  6  a4 

6  -  -  4  +  5.2  a  -  S.2  a3  +  4  a4 

e  “  1  •  2,6  a  +  3.4  a2  •  2.6  a3  +  a4 

It  should  be  noted  that  Equation  (32)  is  the  filter  realized  i.- 
direct  form  and  may  not  be  the  best  to  use  for  the  processing  of  actual 


data.  The  subject  of  how  to  realize  the  filter  as  well  as  more  details 
on  the  s  and  z  transform  are  covered  in  many  texts  and  the  interested 
reader  is  referred  to  some  of  them  [11,12,13,14]. 

Equation  (32)  indicates  that  the  present  filter  output  Y^(nT)  is 
related  to  the  last  four  filter  outputs  as  well  as  to  the  last  four 
measurements  as  well  as  the  present  measurement  Y^fnT) .  It  is  this 
dependence  on  past  measurements  as  well  as  past  outputs  that  introduces 
the  increasing  phase  difference  between  the  input  and  output.  The  order 
of  the  filter  determines  how  many  past  data  values  are  to  be  used.  It 
is  for  this  reason  that  the  order  of  the  filter  should  be  kept  relatively 
small.  Also  the  lower  the  filter  order,  the  easier  it  is  to  program. 

The  results  obtained  when  noisy  data  were  passed  through  the  digital 
filter  just  presented  are  shown  in  Figure  8.  This  figure  shows  the  true 


Figure  8.  Results  of  Processing  Noisy  Data  through 
the  Digital  Filter 
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signal,  «p  the  measured  signal,  and  the  filtered  output,  mf. 

It  should  be  noted  that  the  filtered  output  is  a  large  improvement  over 
the  raw  data  in  indicating  the  true  signal.  The  phase  difference  is 
also  evident  in  this  figure. 

Results  and  Conclusions 

The  results  obtained  when  the  noisy  data  is  first  passed  through 
the  digital  filter  and  then  through  the  adaptive  estimator  are  presented 
in  Figure  9.  In  this  figure  we  see  the  actual  heading  angle,  Ap,  the 


Figure  9.  Results  of  Processing  Noisy  Data 
through  both  Filters 
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actual  measurement,  m^,  the  filtered  measurement  mf  and  the  estimated 
heading  angle  Ap.  Comparing  this  figure  with  Figure  5  it  is  noted  that 
digital  processing  of  the  measurement  is  indeed  an  inprovement  over 
passing  the  raw  data  directly  into  the  adaptive  estimator. 

Using  actual  gyrocompass  data  has  been  investigated  and  some  results 
are  presented  in  Figure  10.  This  figure  shows  the  results  of  passing 
normalized  raw  gyrocompass  data  through  the  digital  filter  as  well  as 


NORMALIZED  TIME  SCALE 


Figure  10.  Results  of  Processing  Actual  Gyrocompass  Data 
through  the  Digital  Filter 

depicting  the  normalized  raw  data.  No  statement  can  be  made  as  to  how  well 
this  approximates  the  true  heading  angle  because  this  fact  is  not  con¬ 
tinuously  known.  The  passing  of  this  digitally  filtered  data  through  the 
adaptive  filter  could  not  be  accomplished  because  the  concurrent  tilt  and 
tilt  rate  i.  'ormation  was  not  available. 
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Based  upon  the  results  obtained  using  the  simulated  data,  where  a 
comparison  can  be  made  to  the  truth,  it  appears  that  the  proposed  process 
is  able  to  satisfactorily  track  a  variable  heading  while  simultaneously 
solving  for  error  model  coefficients.  This  system  can  be  satisfactorily 
used  for  real  time  data  processing  because  the  total  time  required  to 
digitally  filter  and  then  to  adaptively  estimate  the  parameters  is  much 
less  than  the  time  between  sanfles.  Consequently,  this  analytical  tech¬ 
nique  should  be  of  great  value  in  the  determination  of  a  continuous  heading 
reference  that  can  be  used  as  a  laboratory  standard. 

Theoretically,  this  procedure  can  be  modified  to  further  improve  the 
azimuth  measurement  accuracy  of  a  gyrocompass ing  system.  The  provision 
for  the  incorporation  of  additional  measurable  error  model  terms  allows 
for  this  amelioration. 
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